Influence of polydispersity on the critical parameters of an effective potential model 

for asymmetric hard sphere mixtures 



o 
o 

(N 

o 

a 

> 

-i— > 
o3 



G 

O 

o 



> 

On 



O 

O 

-i— > 



c 

o 
o 



X 



Julio Largo 1 and Nigel B. Wilding 1 

department of Physics, University of Bath, Bath BA2 7 'AY, United Kingdom 

(Dated: February 4, 2008) 

We report a Monte Carlo simulation study of the properties of highly asymmetric binary hard 
sphere mixtures. This system is treated within an effective fluid approximation in which the large 
particles interact through a depletion potential (R. Roth et al, Phys. Rev. E62 5360 (2000)) designed 
to capture the effects of a virtual sea of small particles. We generalize this depletion potential to 
include the effects of explicit size dispersity in the large particles and consider the case in which 
the particle diameters are distributed according to a Schulz form having degree of polydispersity 
14%. The resulting alteration (with respect to the monodisperse limit) of the metastable fluid-fluid 
critical point parameters is determined for two values of the ratio of the diameters of the small and 
large particles: q = a a /at = 0.1 and q = 0.05. We find that inclusion of polydispersity moves the 
critical point to lower reservoir volume fractions of the small particles and high volume fractions 
of the large ones. The estimated critical point parameters are found to be in good agreement with 
those predicted by a generalized corresponding states argument which provides a link to the known 
critical adhesion parameter of the adhesive hard sphere model. Finite-size scaling estimates of the 
cluster percolation line in the one phase fluid region indicate that inclusion of polydispersity moves 
the critical point deeper into the percolating regime. This suggests that phase separation is more 
likely to be preempted by dynamical arrest in polydisperse systems. 
PACS numbers: 64.70Fx, 68.35.Rh 



I. INTRODUCTION AND BACKGROUND 

Highly asymmetric mixtures of hard spheres have long 
served as a prototype model for systems of large colloidal 
particles dispersed in a sea of smaller colloids. A key 
physical feature of such systems is the mediation by the 
small particles of so-called depletion forces between the 
large ones [1] . This force has its origin in entropic effects 
associated with the dependence of the free volume of the 
small particles on the degree of clustering of the large 
ones. Although, the typical range of the forces is rather 
limited (of order the diameter of the small particles), they 
can be very strong. 

A longstanding issue in this context concerns the abil- 
ity of depletion forces to engender phase transitions in 
binary hard sphere mixtures. Biben and Hansen [2] ad- 
dressed this matter using integral equation theory, and 
predicted that for sufficient size asymmetry, depletion 
forces engender a fluid-fluid spinodal instability. Other 
theoretical studies have arrived at often conflicting con- 
clusions in this regard (see discussion in ref. [3]), but 
experiments on colloidal systems (eg. ref. [5]), do appar- 
ently confirm a transition, although in certain circum- 
stances it is found to be metastable with respect to a 
broad fluid-solid coexistence region. 

Ideally one would like to settle the matter of the exis- 
tence of a fluid-fluid transition (as well as its stability or 
otherwise with respect to the fluid-solid boundary), by 
computer simulation. Unfortunately, direct simulation 
studies of very asymmetric additive mixtures are severely 
hampered by extremely slow relaxation. Accordingly, all 
such studies to date have been restricted to mixtures of 
relatively low asymmetries, for which a fluid-fluid tran- 
sition is less likely to be observable. While recently de- 



veloped novel algorithms [6, 7] offer some hope of future 
progress in accessing greater size asymmetries, no direct 
evidence has (to date) been obtained for the existence 
of a fluid-fluid phase separation in additive hard sphere 
mixtures. 

In view of these difficulties, a fruitful alternative to 
simulations of the full two component mixture is to at- 
tempt to map it onto an effective one-component sys- 
tem which can be simulated more easily. This "effec- 
tive fluid" approach was taken by Dijkstra, van Roij and 
Evans [3], and by Almarza amd Enciso [4] who proposed 
a model depletion potential by tracing out from the par- 
tition function the degrees of freedom associated with 
the small particles. The resulting interparticle potential 
for the large particles is parameterized by the size ratio 
q = (7 s /<7b between small and large particles, and a reser- 
voir volume fraction of the small particles fj s . Simulation- 
based free energy measurement of the resulting system 
yielded, for values of q < 0.1, a fluid-fluid separation that 
was metastable with respect to a broad solid-fluid coexis- 
tence region. Explicit simulations of the two component 
mixture confirmed the accuracy of the model depletion 
potential as far as the location of solid-fluid and solid- 
solid transitions was concerned, but could not access the 
likely region of fluid-fluid separation. Subsequent work 
by Roth, Evans and Dietrich [8] yielded a more accurate 
depletion potential by fitting to accurate DFT predic- 
tions. However to our knowledge the latter potential has 
to date not been used to study phase behaviour. 

Real colloidal fluids are polydisperse, that is their con- 
stituent particles exhibit an essentially continuous range 
of size, shape or charge. Introducing polydispersity into 
model fluids is known to alter fluid-fluid critical point pa- 
rameters [9, 10] as well as freezing boundaries [11, 12]. It 
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is therefore pertinent to enquire as to the effect of poly- 
dispersity on the location of the metastable fluid-fluid 
transition of hard sphere mixtures. Indeed this question 
has previously been addressed in part by Warren [13] who 
applied the moment free energy method [14] to study a 
polydisperse version of the equation of state of Boublik 
and Mansoori [15] for binary hard sphere mixtures. The 
results showed that for sufficiently large size ratio and 
polydispersity of the large particles, a fluid- fluid spin- 
odal appears in the model. The transition was predicted 
to become more stable with increasing degree of polydis- 
persity. 

Most other theoretical investigations of polydispersity 
in hard sphere mixtures [16-18] have focussed on the form 
of the depletion potential and did not explicitly consider 
the consequences for phase equilibria. The sole study 
of phase behaviour to date (of which we are aware) is 
that of Fasolo and Sollich [11] who applied the moment 
free energy method to the Asakura-Oosawa (AO) model 
[19]. This model describes the limit of maximum non- 
additivity of the small particles and in contrast to addi- 
tive mixtures, the monodisperse AO model is known to 
exhibit a stable fluid-fluid phase transition for size ratios 
q < 0.5 [20, 21]. The introduction of size polydispersity 
to the large spheres [11] was observed to disfavor both 
fluid-fluid and fluid-solid phase separation, though the 
effect was larger for the latter transition. The net result 
was a lowering of the q value necessary for occurrence of 
stable fluid-fluid phase separation, and an increase in the 
stability of this transition with respect to freezing. 

In the present work we apply specialized Monte Carlo 
simulation techniques to an effective fluid model for ad- 
ditive mixtures with a view to elucidating the effect of 
large particle (colloidal) polydispersity on the parameters 
of the fluid-fluid critical point. Our results indicate that 
polydispersity shifts the critical point to lower reservoir 
volume fractions r) s of the smaller particles, and to higher 
volume fractions r/b of the large ones. A determination of 
the percolation threshold using finite-size scaling meth- 
ods shows that the addition of polydispersity moves the 
fluid-fluid critical point deeper into the percolation re- 
gion. We further find that accurate predictions for the 
critical reservior volume fraction f)" lt can be obtained 
by matching the second virial coefficient of the depletion 
potential to that of the adhesive hard sphere model at its 
(independently known) critical point. 



II. MODELS 

A. Depletion potentials 

Two model depletion potentials are considered in this 
work. The first, on which we shall focus primarily, is due 
to Roth, Evans and Dietrich [8], and we shall refer to it 
as the RED potential. It derives from accurate density 
functional theory (DFT) studies of hard sphere mixtures 
and takes the form: 
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Here Rb and R s are the radii of the large and small 
spheres respectively. The scaled depletion potential W = 
W(x, r] s ) is a function of x = h/a s , the distance from con- 
tact measured in units of the small sphere diameter, and 
the volume fraction r/ s of the small spheres.. The param- 
eter e takes the value e = 2 for wall-sphere interactions 
and e = 1 for sphere-sphere interactions. 

Between contact at x = and the location of the first 
maximum x the scaled depletion potential is expressed 
in terms of a cubic polynomial: 



(3W(x,T] s ) = a(ri s ) + b(ri 8 )x + c(ri s )x 2 +d(rj s )x 3 , x < x ■ 

(2) 

The coefficients a,b,c and d were obtained by Roth et al 
by fitting to depletion potentials calculated within DFT. 
For x > x , they assume that the asymptotic regime al- 
ready sets in. For the interaction between two spheres (a 
different expression applies to the asymptotic behaviour 
between a sphere and a hard wall) this is 



pW(x lVs ) = /?W af * m (a;,ffe) ) x>x 



(3) 



where 



l3W asym (x,r) s ) = 



A p (i] s 



^ exp(-a (i] s )<7 s x) (4) 
a s [Ob + h) 

x cos[ai(r] s )a s x - Q p (a s )], x > x 



Here the denominator measures the separation between 
the centers of the spheres in units of a s , and the prefac- 
tors ao(rj s ) and ai(?7 s ) can be calculated from the Percus- 
Yevick bulk pair direct correlation function [8] . The am- 
plitude A p {rj s ) and phase Q p (r] s ) are chosen such that 
the depletion potential and its first derivative are contin- 
uous at xo. They are weakly dependent on the size ratio. 
Fig. 1 shows the form of the potential for the size ratio 
q = 0.1 at a selection of values of rj s . 

The RED potential was tested in ref. [8] by comparing 
with computer simulation results for hard sphere mix- 
tures and was found to perform well for volume fractions 
of the small particles in the range < rj s < 0.3. In 
the simulations to be described below, we use a trun- 
cated version of the RED potential, cutoff at x = 0.3, 
and with no correction. Furthermore we shall employ 
the potential at finite volume fractions rib of the large 
particles, but assume the potential remains two-body in 
form, being parameterised by the reservoir volume frac- 
tion of small particles ry s , which plays a role similar to a 
chemical potential. Clearly f) s — > rj s in the limit 175 — > 0. 
We note that approximate expressions exist which allow 
one to convert from f) s to rj s at finite 775 [3], at least in 
the monodisperse case. 
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FIG. 1: The form of the RED potential for q = 0.1 at a 
selection of values of n a . 



B. Incorporating polydispersity 

The key to incorporating polydispersity into the above 
framework is to generalize the form of the depletion po- 
tential for the case of two interacting large particles of 
different radii R\ and R 2 . This is readily achieved by 
appeal to the Dcrjaguin approximation [24, 25], which 
relates the depletion force between two spheres of differ- 
ent radii (i?i and R2) to the potential between two flat 
plates: 

F ss = 2n RlR l U ww {h) . (6) 
ti\ + il2 

The approximation also yields the force between a sphere 
(of radius R) and a wall: 



The second potential that we have studied, albeit to 
a lesser extent and solely in the monodisperse context, 
is that due to Gotzelmann et al [22, 23]. In contrast to 
the DFT-based RED potential, this was derived purely 
within the framework of the Derjaguin approximation, 
although it too is expressed as a series expansion. We 
shall employ it in the truncated form studied by Dijkstra, 
van Roij and Evans [3] , and refer to as the DRE potential: 



2nRU ww {h) 
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[3A 2 t) s + (9A + 12A 2 )77 s 2 (5) 
-1< A < 



2q 

(36A + 30A 2 )^ 



where A = x — 1. 

Although both potentials have a qualitatively similar 
form at short range, the DRE potential neglects the cor- 
rect damped oscillatory decay at larger particle separa- 
tions. A comparison of the two potentials for size ratio 
q = 0.1 and rj s = 0.3 is shown in fig. 2. 




FIG. 2: Comparison of the RED and DRE potentials for q ■ 
0.1 and fj s = 0.3. 



Clearly, if the two spheres in the first case have equal 
radii (i?i = R2), then 



nRU ww (h), 



(8) 



giving the well known result F sw — 2F SS . 

Returning to the depletion potential of eq. 1, the above 
considerations prompt one to write: 
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where 



R[ = 
R' 2 = 



Ri + R s 

Rs 
R2 + Rs 
Rs 



(9) 

(10) 
(11) 



It is readily verifiable that in the limiting cases R\ — > 00 
and Ri — > R2, one recovers the expression of Roth et al 
[8] (eq. 1) with e = 2 and e = 1 respectively. 

With regard to the effect of this generalization on the 
parameterized form of the potential, we note firstly that 
the quantity x remains unaffected because it is simply 
the distance from contact in units of a s . However, in the 
asymptotic part (eq. 4), the denominator a^ 1 (ai, + h) 
measuring the separation between the centers of the 
spheres in a s units is given in the polydisperse case by 
(jJ 1 [{(Ji +<T2)/2 + /i]. Additionally, while in the monodis- 
perse limit the location of the first maximum of the po- 
tential (at which the asymptotic behaviour is presumed 
to set in) is given by xq = aj l {a s + ct;,), in the polydis- 
perse case one has instead xq = <t~ s ~ 1 [(j s + (<ti + a 2)/ 2]. 
Fig. 3 gives some examples of the influence of dissimilar 
particle sizes on the depletion potential. 

In the present work, we address the situation in which 
the large particles exhibit a continuous variation of sizes. 
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FIG. 3: The form of the RED depletion potential for four 
combinations of the pair radii. In all cases we have set a 3 = 
0.1 
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FIG. 4: The imposed form of f(<7b), corresponding to a Schulz 
distribution (eq. 12) with z = 50. The diameters Ob of the 
large particles are measured in units of cx& = 1 



In order to quantify the form of the polydispersity, we 
label each particle by the value of its diameter <jb- The 
system can then be described in terms of a density distri- 
bution p{(Jb) measuring the number density of particles 
of each o\,. Experimentally, the distribution of colloidal 
particles sizes in a system is (in general) fixed by the 
synthesis of the fluid. To reflect this situation in our 
simulations, we assign p{ob) an ad-hoc prescribed func- 
tional form, which we choose to be of the Schulz type [26] 
defined by the normalized distribution function: 



/(<7b) = T\ (^r) ^ exp 



z + l 

Ob 



■ (12) 



Here z is a parameter which controls the width of the 
distribution, while (Tb = 1 sets the length scale. We have 
elected to study the case z = 50, for which the corre- 
sponding form of /fa) is shown in fig. 4. The associated 
degree of polydispersity is defined as the normalized stan- 
dard deviation of the size distribution: 



(fa 



°b) }'■ 



(13) 



For the Schulz distribution one finds 5 — 1/y/z + 1. With 
z = 50, this formula yields S w 14%. Note that for 
computational convenience, lower and upper cutoffs were 
imposed on the range of allowed particle sizes ov These 
were chosen such that 0.5 < o~b < 1.5. 

The imposed density distribution is related to /(<7b) by 



Pfa) = p£/fa) 



(14) 



where p® is the average number density of large parti- 
cles. Since f(<Jb) is fixed, the form of p(o&) is parame- 
terized solely by p®, variations of which correspond (at a 
given rj s ) to traversing a "dilution line" in the full finite 



dimensional phase diagram [27]. Although this parame- 
terization provides the operational basis for scanning the 
dilution line, we shall (in accordance with convention) 
quote our results in terms of the overall volume fraction 
of the large particles. The latter is related to the density 
distribution via 



f°° 7T 

Vb = J -o-lp{o- b )do- b . (15) 

Finally in this section, we note that within the poly- 
disperse context, the distribution of sizes of the large 
particles implies that the size ratio q can only be defined 
in terms of an average. Accordingly we take q = a s /ab- 



III. COMPUTATIONAL METHODS 

Monte Carlo simulations were performed within the 
grand canonical ensemble GCE using the methods de- 
scribed in refs. [10, 28-30]. Here we briefly outline prin- 
cipal elements of the strategy and refer the interested 
reader to those papers for a fuller description. 

Within the GCE framework, the density distribution 
p(o-b) is obtained as an ensemble average over an instan- 
taneously fluctuating distribution. The form of p{o~b) is 
controlled by the conjugate chemical potential distribu- 
tion p(ab), which was tuned (cf ref. [30]) at all points 
in the phase diagram such as to yield the desired Schulz 
shape f(a) (eqs. 12, and scale p° 14). This tuning was 
achieved by joint use of the non-equilibrium potential re- 
finement (NEPR) method [29], coupled with histogram 
extrapolation [31] in terms of p(ab)- It should be noted, 
however, that the DRE and RED depletion potentials 
do not lend themselves to histogram extrapolation with 
respect to the model parameter i) s which controls the 
form of the interaction potential. This is because f) s does 
not appear as an overall scale factor in the Hamiltonian, 
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a situation which contrasts, for example, to tempera- 
ture reweighting in simpler potentials such as Lennard- 
Jonesium. Consequently in order to scan the phase dia- 
gram with respect to changes in fj s , separate simulations 
were utilized in each instance. 

Our principal aim is a determination of the 
polydispersity-induced shifts in the critical point param- 
eters of the model depletion potentials. To this end we 
have employed a crude version of the finite-size scaling 
(FSS) analysis described in ref. [32]. The analysis in- 
volves scanning the range of p , and fj s until the observed 
probability distribution of the fluctuating instantaneous 
volume fraction of large particles p(f]b), matches the inde- 
pendently known universal fixed point form appropriate 
to the Ising universality class in the FSS limit. Owing 
to the relatively large depth of the interparticle poten- 
tial well (see fig. 2) compared to eg. the Lennard-Jones 
(LJ) potential, the acceptance rate for particle insertions 
and deletions was found to be very low, resulting in ex- 
tended correlation time for the density fluctuations. Con- 
sequently we were able neither to study a wide range of 
system sizes nor obtain data of sufficient statistical qual- 
ity to permit a more sophisticated FSS analysis. Never- 
theless it transpires that our estimated uncertainties on 
the critical point parameters are sufficient to resolve the 
polydispersity-induced trends in the critical point param- 
eters that we set out to identify. 



IV. SIMULATION RESULTS 

Our results are divided into three sections. Firstly we 
locate the fluid-fluid critical point for both the RED and 
DRE potentials in the monodisperse limit. Moving on 
to the polydisperse case, we determine the effect of the 
added polydispersity on the critical point parameters. Fi- 
nally we use finite-size scaling to estimate the locus of 
the cluster percolation threshold in both the mono- and 
poly-disperse cases. 



A. Monodisperse limit 

1. Critical region 

As outlined above, we have tuned the values of fj s 
and /j, until the measured probability distribution of the 
volume fractions of the large particles matched (as far 
as possible given the computational complexity of this 
problem) the universal Ising fixed point form. Fig. 5 
shows distributions obtained in this way for the case of 
the RED potential with q = 0.1 at fj s = 0.3200 and 
f] s = 0.3190. Although the statistical quality is not par- 
ticularly good, comparison of the forms of the distribu- 
tions with that of the fixed point form indicates that 
the given values of fj s straddle criticality, permitting the 
estimate fjf' lt = 0.3195(5). This estimate for the RED 
potential critical point, together for that for q = 0.05, 



RED potential 



'As 



v c b rit 



q 

0.1 0.3195(5) 0.274(10) 
0.05 0.1765(5) 0.289(15) 



DRE potential 



q 

0.1 



V. 



Vb 



0.255(15) 0.286(15) 
[0.289] [0.223] 

0.05 0.151(1) 0.271(15) 
[0.165] [0.235] 

TABLE I: Estimates of critical point parameters, obtained 
using the methods described in the text. Values estimated 
from the data of ref. [3] are given in square brackets. 



and the corresponding estimates for the DRE potential 
are presented in tab. I. 
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FIG. 5: Estimates of the order parameter distribution "p(r\b) at 
fjs = 0.3200 and f) s = 0.3190, for V = (5.2a) 3 . Representative 
error bars are shown. Also included is the fixed point Ising 
magnetisation distribution. All distributions are scaled to 
unit norm and variance via the non-universal scale factor ao- 

With regard to the results, of tab. I, we note that for a 
given potential form, the estimates of r^ rlt appear rather 
insensitive to the value of q. We further note that for a 
given q there is a substantial shift in fj™ 1 * between the 
two forms of the depletion potential. The latter finding 
is perhaps not too surprising given the significant differ- 
ence in the contact value and range of the well depth of 
the RED and DRE potentials, as well as the rather rad- 
ical truncation made by the DRE potential, of the long 
ranged oscillatory part of the interactions (cf. fig. 2). 
Indeed the sensitivity of phase behaviour to the deple- 
tion potential well depth and range has been emphasised 
by Germain et al [33], albeit in the context of fluid-solid 
coexistence. 

We also note significant discrepancies between our es- 
timates of the critical point and those of Dijkstra et al 
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[3] for the DRE potential. Although no error bars are 
quoted in ref. [3] , it seem likely to us that this discrepancy 
is statistically significant. Its source may be traceable to 
the use in [3] of indirect free energy measurements to ob- 
tain the phase diagram, in contrast to the generally more 
accurate direct grand canonical FSS approach employed 
here. Interestingly, our estimates for the critical i) s val- 
ues lie much closer than those of ref. [3] to the results of 
a computation using integral equation theory of both the 
depletion potential and its phase behaviour [34]. 



metastable for a period of time sufficient for us to col- 
lect useful data in the critical region. Unfortunately, the 
freezing became unmanageable when we attempted to ob- 
tain data in the fluid-fluid coexistence region. As noted 
by other authors [3, 35] the coexistence curve of deple- 
tion potentials appears to be rather flat near the critical 
point. Thus even a modest excursion into the two phase 
region results in high liquid densities, which in our expe- 
rience froze very rapidly. 
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FIG. 6: Time evolution of the system volume fraction near 
the parameters of the metastable critical point. Eventually, 
the system freezes spontaneously. 

In ref. [3] it was demonstrated via free energy mea- 
surements that the fluid-fluid critical point for the DRE 
potential is metastable with respect to freezing. While 
we have not attempted to perform a systematic study 
of the freezing transition in the present work, our simu- 
lations confirm the metastability in so far as some runs 
were observed to freeze into an f.c.c. crystal structure. 
An example of the time evolution of the density in such 
a run for the DRE potential is shown in fig. 6. Such 
freezing was also observed for the RED potential indi- 
cating that here too the critical point is metastable with 
respect to crystallization. As an aside, we note that from 
a computational point of view, our observation of freez- 
ing within the grand canonical ensemble is somewhat re- 
markable since the algorithm becomes very inefficient at 
crystal densities. The key factor in achieving this in the 
present case is the inclusion-alongside the standard inser- 
tions, deletions and resizing moves-of particle displace- 
ment moves. Without the latter, the system was not ob- 
served to crystallize on simulation time scales. Another 
apparent factor controlling the ease of freezing appears 
to be whether the crystal lattice parameter is commensu- 
rate with the choice of system box size. We further note 
that our frozen structures do not attain the near-close 
packing densities observed in ref. [3]. This is due to the 
presence of defects in the frozen configurations. 

Notwithstanding the eventual relaxation to a crys- 
talline state, our systems were usually found to remain 



B. Polydisperse case 

Turning now to the polydisperse case, we have ob- 
tained the critical point parameters in a manner simi- 
lar to that employed in the monodisperse limit. Fig. 7 
shows the comparison of our estimates for the critical 
point distribution p(r]b) for the RED potential in both 
the mono and polydisperse cases with q = 0.1. Clearly 
the distribution for the polydisperse case is substantially 
shifted to higher volume fractions compared to that for 
the monodisperse case. Owing to the slow fluctuations of 
rjb, the true average of these distributions could not be de- 
termined to high precision. However, the peak positions 
are rather insensitive to the fluctuations, and on the ba- 
sis of critical point universality one expects that given 
sufficient statistics, the form of the distributions should 
become symmetric. One can therefore estimate r)" lt from 
the average of the peak positions. The results of so doing 
are summarized in tab. II, from which one discovers that 
the principal influence of polydispersity on the critical 
point parameters is a significant decrease in rjg Tlt with 
respect to its monodisperse value, and a significant in- 
crease in ?7g rlt . For the form and degree of polydispersity 
we have studied, the decrease in r)™* is about 6%, while 
the concomitant increase in 77™* is about 17%. 




FIG. 7: Comparison of the distribution p(rit) at the estimated 
critical parameters in the monodisperse and polydisperse sys- 
tems. In both cases the size ratio q = 0.1 and the system size 
is V = (6a) 3 
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RED potential 

a ri crit r> crit 

0.1 0.300(1) 0.336(15) 
0.05 0.1655(5) 0.345(5) 

TABLE II: Estimated critical point parameters for the model 
polydisperse system described in sec. IV B 



The influence of polydispersity on the near critical 
point phase behaviour is further observable in terms of 
particle size fractionation effects. Specifically, when fluc- 
tuations of the instantaneous value volume fraction r)b ex- 
ceed their average value, the distribution of particle sizes 
is shifted to larger diameters; and conversely for fluctua- 
tions of r)b to values lower than the average. The scale of 
the effect is shown in fig. 8 for q — 0.1 at the estimated 
critical point parameters. The presence of such fraction- 
ation implies that the critical point need not lie at the 
apex of the cloud curve that marks the onset of phase sep- 
aration [27]. Indeed we did observe some evidence for a 
weak separation of cloud and shadow curves at rj s = r^"*, 
although precise quantification of the effect was compli- 
cated by a noticibly increased tendency of the polydis- 
perse system to relax to a high density state following a 
fluctuation to high density. The nature of this relaxation 
closely resembled that observed in the monodisperse case 
(cf. fig. 6) and indeed visualization of the arrested con- 
figurations revealed some evidence of crystalline order, 
albeit with a high concentration of defects. We caution, 
however, that these findings should not be interpreted as 
providing strong evidence for a freezing of the polydis- 
perse system because it was not possible to ensure that 
the overall density distribution remained on the dilution 
line during the relaxation to the high density state. 
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FIG. 8: The normalized distribution of particle sizes for 
instantaneous volume fractions r)b below (dashed line) and 
above (dotted line ) the average value, close to the critical 
point. Also shown (solid line) in the overall Schulz "parent" 
distribution 



C. Percolation threshold 

Percolation is a necessary, though not sufficient con- 
dition for gelation and dynamical arrest in colloidal sys- 
tems. Since gelation can affect the ability of experiments 
to observe equilibrium phase behaviour in general, and 
specifically fluid-fluid phase separation, it is important 
to determine the location of the percolation line in the 
phase diagram and its relationship to the fluid-fluid crit- 
ical point. Additionally, it is of interest to ask to what 
extent this relationship is affected by polydispersity 

In order to locate the percolation threshold, it is nec- 
essary to identify pairs of particles that are 'bonded' and 
check for spanning of clusters of such particles. However, 
in contrast to lattice models or fluid systems such as the 
adhesive hard sphere model, the definition of a 'bond' 
in systems with continuous potentials is somewhat am- 
biguous. We therefore adopt a criterion which derives 
from that used for cluster identification in spin models. 
Specifically, we determine the interaction energy u be- 
tween all pairs of particles and assign a bond with prob- 
ability pbond = 1 — cxp(/3u). Clusters of bonded particles 
are then identified using the algorithm of Hoshcn and 
Kopelman [36]. In rcf. [37] Miller and Frenkel identi- 
fied the percolation threshold with those values of the 
model parameters for which the proportion of configu- 
rations containing a spanning cluster is 50%. However, 
finite-size scaling arguments [38] show that better esti- 
mates may be obtained by examining the finite-size be- 
haviour of plots of the fraction of spanning clusters as 
a function of r]b. An example of such a plot is shown 
in fig. 9 for the RED potential in the monodisperse case 
for fj s = 0.28. Data are shown for 4 system sizes, and 
indicate that there is a well defined intersection point at 
77b w 0.21. This intersection point provides a good mea- 
sure of the percolation threshold in the thermodynamic 
limit [38]. By contrast, application of the 50% criterion 
to data for a single system size can considerably overes- 
timate the percolation threshold. 

Percolation lines were determined using this intersec- 
tion method for the monodisperse and polydisperse RED 
potential at q = 0.1. They are shown in fig, 10 together 
with our estimates of the critical point parameters. One 
sees that in both cases the critical point lies well within 
the percolation regime, though much more so for the 
polydisperse system than for the monodisperse system. 



V. LINKING TO THE ADHESIVE HARD 
SPHERE MODEL 

A simple yet general method for finding two potentials 
that are 'equivalent' in a corresponding states sense, is 
to match their second virial coefficient B 2 [37, 39, 40]: 

B 2 = -2ir (V^M - l) r 2 dr . (16) 
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FIG. 9: Fraction of percolating configurations as a function 
of rib for the RED potential in the monodisperse limit. The 
potential parameters are fj s — 0.28, q = 0.1, and data are 
shown for four system sizes. 
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FIG. 10: Percolation line for the RED potential (q = 0.1) for 
both the monodisperse and polydisperse cases, as determined 
by the method described in the text. Statistical errors are 
comparable with the symbol sizes; lines are guides to the eye. 
The system size in both cases was L = (6<r) . Also shown are 
the estimated critical point parameters (cf. sec. II B). 



Here we compare the value of B2 for the RED and DRE 
depletion potentials in the monodisperse limit, with that 
of the adhesive hard sphere (AHS) model [41]. The latter 
comprises hard particles which experience a finite attrac- 
tion only at contact, the strength of which is controlled 
via a 'stickiness' parameter r. The overall interaction 
can be written 



e -/Mr) = e(r _ a) + ° 5{r _ a) t (17) 

lZT 

with r the separation of particle centers and a the particle 
diameter. The second virial coefficient follows as 



qAHS 



3 1 4r 



(18) 





Monodisperse RED potential 




q 


fjs predicted 


fjs simulation 


0.1 


0.320 


0.3195(5) 


0.05 


0.177 


0.1765(5) 




Monodisperse DRE potential 




q 


rj s predicted 


fjs simulation 


0.1 


0.256 


0.255(15) 


0.05 


0.151 


0.151(1) 




Polydisperse RED potential 




q 


fj s predicted 


fjs simulation 


0.1 


0.310 


0.300(1) 


0.05 


0.172 


0.1655(5) 



TABLE III: Comparison of the simulation estimates of the 
critical point parameters with the predictions arising by 
matching the second virial coefficient to that of the critical 
AHS model. 



The AHS model exhibits a fluid-fluid phase transi- 
tion, the critical point of which has been estimated to 
occur [42] at r c = 0.1133(5), p c = 0.508(10). This 
value of r c implies that for the AHS model at criticality, 
B«it = _4.826w with v = (4/3)ttct 3 . It is therefore of 
interest to assess whether, via the matching of B2 values 
for the RED and DRE potentials to that of the critical 
AHS model, reasonable predictions can be made for the 
critical point parameters of the depiction potentials. To 
this end we have numerically evaluated B2 across a range 
of fj s values for each depletion potential and q value of 
interest. By so doing we could determine that value of 
fjs for which B2 matched the value B" lt — — 4.826t>o- 
Table. Ill shows the resulting predictions for 77™* for 
two values of q, together with our simulation estimates. 
Clearly, in each instance, the agreement is remarkable. 

One can attempt to extend the above approach to the 
polydisperse depletion potentials. To do so, we first ob- 
tain the contribution to the second virial coefficient for 



interactions between all pairs of species <7j , Uj ■ The over- 
all coefficient for the mixture can then be approximated 
as a weighted average of pairs [43] , where the weight fac- 
tor is the probability of interaction between a pair of 
particles of size i and size j. In our case, this is given by 
the product of the corresponding values of the normalized 
Schulz distribution (eq.12): 



B 2 = 



f(<7 l )f(crj)B2(a l ,(Tj)d<T l daj (19) 



Matching to B"^ — — 4.826^o as before, one obtains 
for the two q values studied, the predictions for ^ rlt 
shown in tab III. Here the agreement with the simula- 
tion estimates is less impressive than in the monodisperse 
case. Although the absolute value of the prediction still 
agrees to within about 3% with the simulation estimate, 
and the sign of the polydispersity-induced shift in 77™* is 
correctly predicted, its magnitude is underestimated by 
a factor of two. 
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The larger relative discrepancy between the predicted 
and measured r)™' may point to a breakdown in the pres- 
ence of polydispersity of the assumed model invariance 
of the critical B 2 value. Indeed one might expect such a 
failure because the value of B 2 is based solely on the pair 
potential and takes no account of the ability of a poly- 
disperse fluid to exploit local size segregation in order to 
pack more effectively than a corresponding monodisperse 
one. In order to address this issue directly, one would 
require estimates of critical point B 2 values for a poly- 
disperse version of the AHS model. To our knowledge no 
simulation estimates of the liquid-gas transition currently 
exist for a polydispcrse AHS model. Indeed, the matter 
is complicated by the fact that there is no unique model 
for polydispersity in such a system. Very recently, how- 
ever, a number of physically reasonable models for poly- 
dispersity in AHS system have been proposed by Fantoni 
et al [44], who investigated the corresponding phase be- 
haviour using integral equation theory. From ref. [44], 
one can deduce that the presence of polydispersity sig- 
nificantly decreases the magnitude of B 2 at the critical 
point compared to the monodisperse limit. This trend in 
B 2 is of the correct sign and overall magnitude to push 
the predictions for rj" lt for our depletion potentials closer 
to the simulation estimates. Unfortunately since no data 
were reported for exactly the same degree of polydisper- 
sity (S — 14%) studied in the present work, no direct 
comparison of B 2 values is possible. 

In an attempt to throw additional light (albeit indi- 
rectly) on the discrepancy between the measured and 
predicted critical point parameters, we have studied the 
effect of introducing polydispersity on the critical point 
B 2 value for the Lcnnard- Jones fluid, which is a com- 
putationally more tractable system than the AHS model 
[45]. The corresponding potential is 



with tij = (Ji<Tje, Uij = {o l + <Jj)/2 and r io = |rj - r,-|. 
The potential was cutoff for > 2.5aij and no tail cor- 
rections were applied. For the monodisperse limit, the 
critical temperature occurs at T c = 1.1876(3) [32] and 
one finds B 2 = — 6.621v , which lies within the range of 
'typical' critical point B 2 values found in surveys of a 
wide range of model potentials [39, 40]. If, on the other 
hand, a is distributed according to a Schulz form (eq.12) 
with z = 50, as used elsewhere in this work, simulations 
yield a critical temperature of T c — 1.384(1), for which 
(using eq. 19), one finds B 2 = — 5.759«o, which is sig- 
nificantly smaller in magnitude that the monodisperse 
value. 

If one now makes the (not unreasonable) assumption 
that given the same form and degree of polydispersity, a 
comparable fractional change of B 2 will ensue in other in- 
teraction potentials, one can estimate the expected crit- 
ical point value of B 2 for the polydisperse AHS model 



(and hence also the polydisperse depletion potentials), 
as B 2 = -4.826 x 5.759/6.621 w -4.2v . The resulting 
predictions for the critical reservoir volume fraction of 
the small particles are fy™* = 0.304, for q = 0.1, and 
ry" 1 ' = 0.169 for q = 0.05, both of which agree within 
error with our simulation results. Of course in view of 
our assumptions, this accord may be fortuitous, but it is 
nevertheless suggestive that a more detailed assessment 
of the effect of polydispersity on critical point B 2 values 
in other systems (specifically the AHS model) would be 
a worthwhile avenue for further study. 



VI. CONCLUSIONS AND DISCUSSION 

To summarize, we have determined the effect of intro- 
ducing polydispersity on the fluid-fluid critical point pa- 
rameters of a model depletion potential for highly asym- 
metric additive hard sphere mixtures. For the particular 
realization of the polydispersity considered, the critical 
point is found to shift (with respect to the monodisperse 
limit) to smaller values of the reservoir volume fraction of 
the small particles r) s and to larger values of the volume 
fraction of the large particle r\\,- It seems reasonable to 
assume that the direction of the shifts should be a general 
trend, common to other functional forms and degree of 
polydispersity. Indeed the same trend has also recently 
been observed in a study of colloidal polydispersity in the 
AO model [11]. 

Beyond this, our results show that inclusion of polydis- 
persity pushes the whole fluid-fluid binodal deeper into 
the percolating regime. Since colloidal fluids are known 
to form a gel [46, 47] for sufficiently high r) s , it would 
seem that the presence of polydispersity increases the 
likelihood that direct observations of fluid-fluid phase co- 
existence is complicated by dynamical arrest. 

Additionally we demonstrated that excellent predic- 
tions for the value of fj$ Tlt follow from matching the sec- 
ond virial coefficient B 2 of depletion potentials to the 
critical B 2 value of the adhesive hard sphere model. The 
quantitative accuracy of the predictions is undoubtedly 
due in large part to the very short ranged nature of the 
depletion potentials; similar studies comparing critical 
point B 2 values for a range of other potentials [39, 40] 
did not find such a high degree of accuracy. Nevertheless 
our observation should prove generally useful in reduc- 
ing the effort required to locate criticality in depiction 
potentials. It is intriguing however, that the accuracy of 
the predictions was reduced on incorporating polydisper- 
sity, suggesting that (perhaps due to changes in packing 
ability due to local size segregation effects), the inclusion 
of polydispersity in a model does not leave B 2 invari- 
ant at the critical point. Comparisons of the critical B 2 
value for a monodisperse and polydisperse LJ fluid con- 
firmed a significant difference in this regard. Moreover, 
the magnitude of the effect was sufficient to account for 
the discrepancy in the observed and predicted 77™* for 
the polydisperse depletion potential. Clearly, however, 
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further work is called for in order to elucidate this mat- 
ter more fully. 

Obviously knowledge of the shift in the critical point 
parameters is in itself insufficient to determine whether 
polydispersity renders the fluid-fluid transition stable 
with respect to fluid-solid coexistence. Although we did 
observe a spontaneous relaxation of the near critical poly- 
disperse system to a high density state showing some 
crystalline order, this finding should be treated with cau- 
tion because the system departs from the dilution line 
during the formation of the new state. It would thus be 
interesting in future work to try to study explicitly the 
effects of polydispersity on the freezing transition. As 
well as providing assessment of the overall stability of the 
fluid-fluid critical point, freezing in polydisperse fluids is 
a matter of considerable interest in its own right. Indeed 
recent theoretical calculations for the AO model indicate 
an increasing richness of fluid-solid and solid-solid phase 
behaviour as the degree of polydispersity is increased [11]. 
To tackle this computationally, however, is a considerable 
challenge, but one which might be met by extending to 
polydisperse system novel computational methods which 
have hitherto only be deployed in the monodisperse con- 



text [48, 49]. 

Finally, we remark that while the present work has con- 
sidered solely the case of polydispersity of the large par- 
ticles, the converse situation of small particles polydis- 
persity is clearly of interest and practical relevance too. 
While there have been several theoretical studies which 
have considered this scenario (see eg. refs. [17, 50]), we 
know of no simulation studies to date. An extension 
of the methods utilised here to address this case would 
doubtless be a worthwhile endeavor-one which we hope 
to undertake in future work. 
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